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ABSTRACT 

The terrestrial fossil record shows a significant variation in the extinction and origination rates of 
species during the past half billion years. Numerous studies have claimed an association between 
this variation and the motion of the Sun around the Galaxy, invoking the modulation of cosmic rays, 
gamma rays and comet impact frequency as a cause of this biodiversity variation. However, some of 
these studies exhibit methodological problems, or were based on coarse assumptions (such as a strict 
periodicity of the solar orbit). Here we investigate this link in more detail, using a model of the Galaxy 
to reconstruct the solar orbit and thus a predictive model of the temporal variation of the extinction 
rate due to astronomical mechanisms. We compare these predictions as well as those of various 
reference models with paleontological data. Our approach involves Bayesian model comparison, which 
takes into account the uncertainties in the paleontological data as well as the distribution of solar 
orbits consistent with the uncertainties in the astronomical data. We find that various versions of 
the orbital model are not favored beyond simpler reference models. In particular, the distribution 
of mass extinction events can be explained just as well by a uniform random distribution as by any 
other model tested. Although our negative results on the orbital model are robust to changes in the 
Galaxy model, the Sun's coordinates and the errors in the data, we also find that it would be very 
difficult to positively identify the orbital model even if it were the true one. (In contrast, we do find 
evidence against simpler periodic models.) Thus while we cannot rule out there being some connection 
between solar motion and biodiversity variations on the Earth, we conclude that it is difficult to give 
convincing positive conclusions of such a connection using current data. 

Keywords: astrobiology — Earth — Galaxy: kinematics and dynamics — methods: statistical 



1. INTRODUCTION 
1.1. Background 

Over the course of Earth's history, evolution has pro- 
duced a wide variety of life. This is particularly apparent 
from around 550 Myr ago — the start of the Phanero- 
zoic eon — when hard-shelled animals first appeared and 
were preserved in the fossil record. Since then we ob- 
serve a general increase in the dive rsity of life, but with 
significant variation superimposed (ISepkoski et al.l 120021 : 
IRohde fc MiHieil [200l lAlrov et al.ll2008D . The largest 
and most rapid decreases in biodiversity — defined here 
are the number of genera extant at any one time — are 
referred to as mass extinctions. 

The cause of these variations in biodiversity in general, 
and mass extinctions in particular, has been the subject 
of intense study and speculation for over a century. Many 
mechanisms have been proposed for the observed varia- 
tion, which we can place into four groups. 

First, the variations are the result of inter-species in- 
teractions. Species compete for limited resources, and 
as one species evolves to compete in this struggle for 
survival, so other species will evolve too. This idea 
has been r e ferred to as th e "Re d Queen hypothesis" 
(jVan ValenI ()1973D : iBentonI (|2009D : in reference to the 
Red Queen's race in Lewis Carroll's Through the Looking 
Glass, where Alice must run just to keep still). Recent 
ecological studies indicate that the interaction between 
species can minimize co mpetition and enhance biodiver- 
sity (Geor ge &: Had '2009). although it is not obvious that 
these biotic factors a re the main cause of large- scale pat- 
terns of biodiversity ()Bentonll2009t IAhovll2008l) . 



Second, the environment changes with time, and 
species will evolve in response to this. This ideas is some- 
times called the "C ourt Jester hypothesis" (Barnosk^ 
'20011; IBentonI 12009( 1. Some of these (abiotic) geo- 
logical changes are relatively slow, such as plate tec- 
tonics, atniospheric composition, and global climate 
(ISigurdssonlll98 8': 'Crow lev fc NorthllT988t IWignal]|l2001l: 
iMartf fc Erns? r2005: Fe ulneri 120091 : iWignall et al.ll2009| l . 
Others may be more rapid. Large-scale volcanism, for 
example, would inject dust, sulfate aerosols and car- 
bon dioxide into the atmosphere, resulting in a short- 
term global cooling, reduced photosynthesis, long peri- 
ods of acid rain, and resulting ultimately in a long-term 

fibal warming (on a timescale of 10^ vr: lMarti fc ErnstI 
)05)). 

Third, extraterrestrial mechanisms could be involved, 
either through a direct impact on life or by changing 
the terrestrial climate. Variations in Earth's orbit 
(mostly its eccentricity) over ten to one hundred 
thousa nd year time-scal e s are responsible for the ic e 
ages (IHavs et all [19761 : IMuller fc MacDonaldl 1200(11) . 
Extraterrestrial mechanisms on longer time scales 
could also pla y a role. These i n clude solar variability 
(IShavivl l200a ILockwoodl l2005l : iLockwood fc FrohlichI 
2007D. asteroid or comet imp acts (|Shoemakeij 119831: 
Alvarez et all [T 9M 'gIm? 1991) cosmic rays (jShavivl 
20051: ISloanfc Wolfcndalc 200|), s upernovae (SNe) 
and ga mma-ray burst (GRBs; [Ellis fc Schramml 
([T995h: iMelott Adrian L. fc Thomas Briaii C.I i20M- 
iDomainko et al.1 ()2013[ ): for a review see iBailer- Jones! 
(2009D). For example, cosmic rays might influence 
Earth's climate if they play a significant role in cloud 



2 



format ion (through the formation of cloud condensation 
nuclei; iCarslaw et all (|2002D : iKirkEvI (|2007D V Secondary 
muons resulting from cosmic rays - as well as high energy 
gamma rays from SNe - could kill organisms directly 
or damage their DNA (iThorsettl llQQSt iScalo fc Wheeled 
[2nnl lAtri fc Melottl[20m 

Finally, the apparent variation may be, in par t, the 
result of uneven preserv ation and sampling bias ()Raud 
[1973 1 Alrovl[l996l: lPetersl[20Ql . Fossilization is relatively 
rare and some animals are more likely to be preserved 
than others. Furthermore, the degree of preservation of 
marine spieces (more common in the fossil record) de- 
pends on the amount of continental outcrop a vailable 
at any time, and th is depends on the sea level ([Hallaml 
[1981 IHollandl[20Tl) . The number of species or genera 
living at any one time is not observed but must be re- 
constructed from the times at which species appear and 
disappear, which implies some kind of sampling or model- 
ing. This can introduce a bias, although it m ay be dimin- 
ished to soine deg ree by various techniques (jAlrov et al.l 
[2MI: [MroA^[20ia) . 

These four types of cause of biodiversity variation are 
not mutually exclusive. They probably all acted at some 
point, and will also have interacted. For example, an as- 
teroid impact could release so much carbon dioxide that 
long-term global warming has the biggest impact on bio- 
diversity. Alternatively, a cool period would lower sea 
levels, leaving less continental shelf for the preservation 
of marine fossils, even though biodiversity itself may be 
unchanged. 

Given the limited geological record, untangling the rel- 
evance of these different causes in most individual cases 
is difficult, if not impossible. More promising might be 
an attempt to identify the overall, long-term significance 
of these potential causes. The goal of this article is to do 
that for extraterrestrial phenomena. 

Many of the astronomical mechanisms mentioned 
above are ultimately caused by the presence of nearby 
stars. Stars turn SNe, the source of gamma rays, 
and their remnants a re a major source of cosmic rays 
(jKovama et al.l [l995| ). Stars perturb the Oort cloud, 
the main source of comets in the inner solar system 
(jRampino fc Stotherd[T98l iGarcia-Sanchez et al.ll2001[ ). 
Broadly speaking, when the Sun is in regions of higher 
stellar density, it is more exposed to extraterrestrial 
mechanisms of biodiversity change. In its orbit around 
the Galaxy (once every 200-250 Myr or so), the Sun's 
environment changes. For example, it oscillates about 
the Galactic plane with a (quasi) period of 50-75 Myr 
([Bahcall &: BahcaIilll985D , and in doing so moves through 
regions of more intense star formation activity in the 
Galactic pla ne. This is particularly true if the Sun crosse s 
spiral arms (|Gies fc Hel sel"2005': 'Leitch fc Vasishtlll998l ). 
which it may do every 100-200 Myr or so. 

Such changes in the solar environment have been 
used as the basis for many claims of a causal 
connection between the solar motion and mass 
extinctions and/or climate change. Typically, au- 
thors have identified a periodicity in the fossil 
record and then connected this to a plausible pe- 
riodic it y in the solar motion (lAlvarez fc Mulled 
T98l IRaup fc SepkosB [l98l: IDavis et al.l 1 1984 : 
Mullei' '1988'; 'Shaviv '2003'; 'Rohde & MiHIei 2005 : 
Mclott Adrian L. fc Banibach Richard k] I2OII : 



iMelott et al] [201I) . These comparisons are fraught 
with problems, however, some of which remained 
unmentioned by the authors. The first is the fact 
that the solar motion and past environment are poorly 
constrained by the astronomical data, s o a wide range 
of p l ausible periods are permis sible (jOverholt et al.l 
120091: iMishurov fc Acharoval 1201 If) , yet the coincident 
one is naturally chosen. Second comes the fact that 
the solar motion is not strictly periodic even under the 
best assumptions. Third, many of these studies have 
not performed a careful model comparison. Typically 
they identify the best-fitting period assuming the 
periodic model to be true, but fail to accept that 
a non -periodic model migh t explain the data eve n 
better dKitchell fc Penal [1981 IStider fc WameJ [1981 . 
In some cases a significance test is introduced to exclude 
a specific noise model, but this is often misinterpreted, 
and the resul ting significance ove restimated. The reader 
is referred to iBailer- Jones! (|2009[ ) for an in-depth review 
and references. 

1.2. Overview 

Here we attempt a more systematic assessment of the 
possible role of the solar orbit in modulating extrater- 
restrial extinction mechanisms. Our approach is new 
in a number of respects, because we: (1) do a numer- 
ical reconstruction of the solar orbit (rather than just 
assuming it to be periodic); (2) take into account the ob- 
servational uncertainties in that reconstruction; (3) use 
models which predict the variation of probability of ex- 
tinction with time (rather than assuming that extinc- 
tion events occur deterministically, for instance); (4) do 
proper model comparison (rather than using p-values in 
an over-simplified significance test); (5) compare not only 
the orbital model with the fossil record but also numerous 
reference (non-orbital) models, such as periodic, quasi- 
periodic, trend, periodic with constant background, etc. 

Our method is as follows. Adopting a model for the 
distribution of mass in the Galaxy, we reconstruct the so- 
lar orbit over the past 550 Myr by integrating the Sun's 
trajectory back in time from the current phase space co- 
ordinates (position and velocity). This gives us a time 
series of how the stellar density in the vicinity of the Sun 
has varied over the Phanerozoic. We then assume that 
this density is approximately proportional to the terres- 
trial extinction probability (per unit time). That is, we 
adopt a non-specific kill mechanism linking the solar mo- 
tion to terrestrial biodiveristy. This is naturally a strong 
and rather simple assumption, but it should be empha- 
sized that we are interested in the overall plausibility of 
extraterrestrial phenomena rather than trying to identify 
a specific cause of individual extinction events. The re- 
sulting time series is then compared to several different 
reconstructions of the biodiversity record. 

A significant source of uncertainty in the reconstructed 
solar orbit is the current phase space coordinates (or "ini- 
tial" conditions). We therefore sample over these to build 
up a set of (thousands of) possible solar orbits, and com- 
pare each of these with the data. The comparison is 
done by calculating the likelihood of the data for each 
orbit. Rather than finding the single most likely orbit, 
we calculate the average likelihood over all orbits. This is 
important, because it properly takes into account the un- 
certainties (whereas selecting the single most likely orbit 
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would ignore them entirely). Indeed, these initial con- 
ditions can be considered as the six parameters of this 
orbital model (for a fixed Galactic mass distribution). 
We are, therefore, averaging the likelihood for this model 
over the prior plausibility of each of its parameters. This 
average — or marginal — likelihood is often called the 
"evidence" . This is just the standard, Bayesian approach 
to model assessme nt, which avoids the various flaws of 
hypo thesis testing peffrevsll 19611: lWinkleiill972l: iMacKavl 
120031 ) , but unfortunately it has seen little use in this field 
of research. 

The next step is to compare this evidence with that 
calculated for various reference models (sometimes called 
"noise" or "background" models, depending on the con- 
text). One such model is a purely sinusoidal model, 
parameterized by an amplitude, period and phase. We 
generate a large number of realizations of the model for 
different combinations of the parameters, calculate the 
likelihood of the data for each, and average the results. 
This averaging plays the crucial role of accommodating 
the complexity of the model. A complex model with lots 
of parameters can often be made to fit an arbitrary data 
set well. That is, it will give a high maximum likelihood. 
But this does not make it a good model, precisely be- 
cause we know that it could have been made to fit any 
data set well! Such models are highly tuned, so while the 
maximum likelihood fit may be very good, a small per- 
turbation of the parameters results in poor predictions. 
Unless supported by the data very well, such models are 
less plausible. A simpler model, in contrast, may not 
give such an optimal fit, but it is typically more robust to 
small perturbations of the model parameters or the data, 
so gives good fit over a wider portion of the parameter 
space. The model evidence embodies and quantifies this 
trade off, which is why it — rather than the maximum 
likelihood — should be used to compare models. 

We have selected four data sets for our study. The first 
two are compil ations of the vari a tion o f ext inction rate 
over t ime, from lRohde fc MulleU ()2005D and lAlrov et al.l 
(|2008f l. In the latter we use the extinction rate stan- 
dardized to remove the sampling bias. Both report a 
magnitude as a function of time. The second two data 
sets just record the time of mass extinction events. Here 
we take the times of the "big 5" mass extinctio ns and 18 
mass extinctions identifi ed bv iBambachI (l2006f) b ased on 
Sepkoski's earlier work (Sepkos ki fc Raupl 119861 ). Each 
mass extinction is represented as a (normalized) Gaus- 
sian on the time axis, the mean representing the best 
estimate of the date of the event and the standard devia- 
tion the uncertainty. We refer to these as "discrete" data 
sets, as they just list the discrete dates at which events 
occur (we do not use any magnitude information). The 
two rate data sets we therefore refer to as "continuous" 
(even though in practice the rates are also recorded at 
discrete time points). 

The paper is arranged as follows. We first introduce 
the data sets. In Section [3] we explain the modeling ap- 
proach, how we calculate the likelihoods and evidence, 
and how we compare the models. In Section |4] we de- 
scribe how we reconstructed the solar orbit, and also 
quantify the degree of periodicity typically present (as 
a strict periodicity has often been assumed in the past). 
In Section [5] we calculate and compare the evidences for 
the various models and data sets and test the sensitivity 



of the results to the model parameters and uncertainties 
in the data. We conclude in Section |6l 

2. PALEONTOLOGICAL DATA 

We adopt four data sets: two discrete time series (se- 
quence of time points with age uncertainties) giving the 
dates of mass extinctions, and two continuous time series 
giving the smoothed and normalized extinction rate as a 
function of time. 

2.1. Discrete data sets 

Sepkoski and others have i dentified five extinctio n 
events to be "mass extinctions" (jSepkoski fc R,auDl[T986[ ). 
often referred to as the "big five". Other studies have 
identified different candidates for these, or have identified 
a "big N" for some other value of N. For example, Bam- 
bach et al. identify t hree mass extinction events as be - 
ing globally distinct ()Bambach Richard K. et ahl 120041) . 
Here we ad opt a set of 18 ma ss extinction events (or B18) 
selected by 'Bambach (2006') using an updated Sepkoski 
genus-level database. They are consistently identifiable 
in different biodiversity data sets and when using differ- 
ent tabulation methods. The second of our discrete data 
sets is the "big five" as identified from among the B18. 
Other choices of events are of course possible and our 
results will, in general, depend on this choice (although 
as we will see the results are rather consistent). 

The times and durations of the events are listed in Ta- 
ble[T] The time, r, is the mid-point between the start age 
and the end age of the substage in which the extinction 
occurred, and the substage duration, d, is the difference 
between these. The geological record does not resolve the 
extinction event, so the extinction presumably took place 
more rapidly than this substage duration. In that case r 
is our best estimate of the true (but unknown) time, t, at 
which the extinction occurred, and c? is a measure of our 
uncertainty in this estimate. Uncertainty is represented 
by probability, so we interpret an "event" as the proba- 
bility distribution P{T\t,d). This is the probability that 
we would measure the event time as r, given t and d. We 
could represent this as a rectangular ("top-hat") distri- 
bution of mean t and width d, but this assigns exactly 
zero probability outside the substage duration, which im- 
plies certainty of the start and end ages. Even though 
their (relative) ages have uncertainties which are less the 
event duration, we nonetheless accommodate some un- 
certainty in these start and end ages by considering each 
event to be a Gaussian distribution with mean t, and 
standard deviation a equal to the standard deviation of 
the rectangular distribution, which is cr = dj \fVl. This 
Gaussian distribution is broader than the corresponding 
rectangular distribution. Note that the intensity of the 
mass extinction is not taken into account. A Gaussian is 
normalized, so its peak value is determined by its stan- 
dard deviation (Figure [TJQ 

2.2. Continuous data sets 

^ One might think that a more natural interpretation of an event 
is P{t\T,d), the probabihty that the true event occurs at t given 
the measurements. But here we are considering the measurement 
model (or noise model), that is, given some true time of the event, 
what possible times might we measure, the discrepancy arising on 
account of the finite precision of our measurement process. 
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Figure 1. Four data sets used in this study. The top row shows the discrete data sets: the B5 (left) and B18 (right) mass extinction 
events. These can be interpreted as the extinction probability density function (PDF), which is proportional to the extinction fraction per 
unit time (i.e. a rate). The bottom row shows the continuous data sets, which give the extinction rate: RM (left) and AOS (right). 



Table 1 

The B18 mass extinction events, with the B5 shown in bold. 
= before present. 

Time (t) / Myr BP Substagc duration (d) / Myr 
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3.53 


489.350 


2.10 
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2.15 


499.950 


2.10 


519.875 


2.75 



The discrete data sets are naturally biased in that they 
BP only select periods of high extinction rate. It may well 
be that extraterrestrial phenomena are only relevant in 
causing (or contributing to) mass extinctions, but a priori 
it is natural to ask how the overall extinction rate varies. 
The extinction rate, E{t), is the fraction of genera which 
go extinct in a stratigraphic substage divided by its du- 
ration. This is directly proportional to the variation of 
extinction probability per unit time. For one of our ex- 
tinction rate data sets we use t he linearized and inter- 
polated data s et construc t ed by iRohde &: Mulled ()2005f ) 
as reported in iBambachI ()2OO60 . We denote this RM. 
The other data set is the "thre e-timer" e xtinction rate 
from the Paleobiology databas60 (|Alrov et a l. 2008). The 
data are binned into 48 intervals averaging 11 Myr in du- 
ration. The counts are derived from 281491 occurrences 
of 18 541 genera within 42 627 fossil collections. We use 
the data set processed using their subsampling method in 
order to reduce the sampling bias, and denote this AOS. 
Both data sets are reported as lists of extinction rates at 
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specific times, {rj,Tj}. These two continuous data sets 
are plotted in the lower row of Figure [TJ 

3. DATA MODELING METHOD 

3.1. Overview of Bayesian model comparison 

The goal of this work is to compare how well vari- 
ous models predict the paleontological data sets. We do 
this by calculating, for a given data set, the Bayesian 
evidence for each model. If the models are equally prob- 
able a priori, then the one with the highest evidence is 
the best predictor of the data. This does not exclude 
the possibility that there exists a better model which we 
have not yet tested. But it at least allows us to conclude 
that the lower evidence models are neither appropriate 
nor sufficient explanations of the phenomenon. (Indeed, 
we never assume a model is "true" , just better than the 
alternatives.) 

The mode l hng a pproach is described in full by 
iBailer- Jones! ()2011altf ). so will only be outlined here. 
Let D denote the paleontolgical time series, and M the 
model. Examples of M are the orbital model or the pe- 
riodic model. Each model has a set of parameters, 9. 
The likelihood of the model, P{D\6, M), is the probabil- 
ity of obtaining D from model M with its parameters set 
to some specific values of 9. Normally we do not know 
the exact values of these parameters, and the data — 
being noisy and imperfectly fit by the model — do not 
determine them exactly either. We therefore average the 
likelihood over all possible values of 6, weighting each by 
how plausible that value of 61 is. This weighted average is 
the evidence. This weighting is given by the prior prob- 
ability distribution, P{9\M). In the case of the orbital 
model, where 9 is the current phase space coordinates of 
the Sun, P{9\M) is determined by the observational un- 
certainties in these measurements. Mathematically the 
evidence is 

P{D\M)^ P{D\9,M)P{9\M)d9 . (1) 
Je 

It gives the probability of getting the data from the 
model, regardless of the specific values of the parame- 
ters, i.e. it measures how well the model explains the 
data. The absolute value of the evidence is not of inter- 
est, so we generally deal with the ratio of two evidences 
for two models, known as the Bayes factor. The evidence 
is a far better measure of the suitability of a model than is 
the p - value jje ffrcvs 1961; Kass & Raftcrvl fTOOSl iWinklei! 
[I97a lMacKavfl2003; .Bailer- Jones,.200a) . 

It is worth stressing again that, by averaging over the 
parameters, the evidence is not sensitive to model com- 
plexity per se. This is in contrast to the likelihood at 
the best fitting parameters (the maximum likelihood): 
a more complex (flexible) model will always fit the data 
better, and so will always deliver a higher maximum like- 
lihood. The evidence reports the average likelihood, so 
it will only increase if the extra complexity gives a net 
benefit over the plausible parameter space. The model 
complexity does not then need to be considered sepa- 
rately in some ad hoc way. 

3.2. The likelihood 
3.2.1. Discrete data sets 



As explained in Section 12. 1[ each event in a discrete 
data set can be considered as a Gaussian distribution. 
For event j this gives the probability that the mass ex- 
tinction occurs at time tj, given that the true time is tj 
and the uncertainty in our measurement is aj. That is. 



In order to compare the measurement of this event with 
the predictions of a model we calculate the likelihood. 
This is given by an integral over the unknown true time 

P{Tj\aj,9,M)= / P{Tj\aj,tj,9,M)P{tj\aj,9,M)dtj 




P{T,\a,,t,)P{t,\9,M)dt^ . (3) 



The first term in the integral - Eqn. [2] - is sometimes 
called the measurement model. The second term is the 
prediction of the time series model, i.e. the probability 
(per unit time) that a mass extinction occurs at time tj . 
Our time series models are, therefore, stochastic in the 
sense that they do not attempt to predict when mass 
extinctions occurred, but rather how the probability of 
occurrence of a mass extinction varies over time. Note 
that the likelihood is just measuring the degree of overlap 
between the data and the model predictions, averaged 
over all time. 

Eqn.[2]gives the likelihood for a single event. Assuming 
all events are measured independentl}0, then the likeli- 
hood for all the data is just the product of the event 
likelihoods 

PiD\9,M) = l[P{T,\<j,,9,M) (4) 

j 

where D = {tj}. For the sake of the likelihood and 
evidence calculation we do not consider the {cTj} as data, 
although they are of course measured. That is because 
D is defined as just those quantities which are predicted 
by the measurement model. 

3.2.2. Continuous data sets 

Mathematically the likelihood for the continuous data 
is very similar as in the discrete case, but the interpreta- 
tion is different. 

Consider the measurement model, P(r|cr, i), in Eqn. [2] 
(we consider just one event so drop the subscript j). We 
have interpreted this as the probability of a discrete ex- 
tinction event being measured at r, but we could equiv- 
alently interpret it as the probability density (i.e. proba- 
bility per unit time) of extinction at time r. Now, rather 
than characterizing the probability density as a Gaussian 
with mean t and standard deviation cr, we could consider 
an arbitrary function, characterized by a series of top-hat 
functions, {pi} (a histogram), each top- hat characterized 
by a center ti, height r^, and width Si. We can then re- 

^ More precisely, the events are assumed independent given the 
model and its parameters. This is probably a reasonable assump- 
tion given that the events are distributed quite sparsely over the 
Phanerozoic, and that the separations between them are generally 
much longer than their substage durations. 
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place P(t|(T, i) with J2iPii'''\^ij^i) where 

, I Jr^when - 6,/2 < t < U + 6^/2 , . 

P^i^l^^^^^) ~ \ otherwise . 

We now see that 

lim Vp.(T|i„5.)-i?W (6) 

i.e. we get a continuous function of the variation of the 
extinction probabihty with time t (as r and t become 
equivalent in this limit). The likelihood is then 

PiD\9, M) = j^E{t)P{t\9, M)Ai . (7) 

In practice we characterize E{t) using the extinction rate, 
rj, tabulated at each time Tj, which is equivalent to 
assuming that extinction rate is constant over the sub- 
stage (or that we have zero uncertainties on the measured 
times) . 

We can actually apply this interpretation to the dis- 
crete data sets too. In both cases, the data provide the 
variation of extinction probability (per unit time) as a 
function of time, or something proportional to that. The 
proportionality constant is irrelevant, because we keep 
the data fixed when comparing different models using 
the evidence. 

3.3. Time series models 

The time series models appear in the equation for the 
likelihood in the form P{t\6, M), i.e. the extinction prob- 
ability (per unit time) as a function of time predicted by 
model M at parameters 0. It is important to realize that 
this probability density function (PDF) over t is normal- 
ized, i.e. integrates over all time to unity. This is key 
to model comparison, because a model which assigns a 
lot of probability to extinctions at some particular time 
must necessarily assign lower probability elsewhere. This 
follows because we are not trying to model the absolute 
value of the extinction rate, but just its relative varia- 
tions. 

In addition to specifying the functional form of the 
models we must also specify the prior probability distri- 
bution of the model parameters, P{9\M). This describes 
our prior knowledge of the relative probability of differ- 
ent parameter settings. For example, given the time scale 
in the data, we are not interested in models with time 
scales less than a few million years or more than a few 
hundred million years. It is often difficult to be precise 
about priors, and the evidence and therefore Bayes fac- 
tors often depend on the choice. The choice of prior must 
therefore be considered part of the model (e.g. "periodic 
model with permissible periods between 10 and 100 Myr" 
is distinct from "periodic model with permissible periods 
between 50 and 60 Myr"). We investigate the sensitivity 
of the results to changes in the prior in Section 15.31 Ex- 
cept for the orbital model, we adopt a uniform prior over 
all model parameters over the range specified in Table [3l 

Table [2] summarizes the functional form of the models, 
which are now briefly described. Figure [2] plots examples 
of some of these models. The range of the data is taken 
to be 0-550 Myr BP. 

Uniform: Constant extinction PDF over the range of 
the data. This has no parameters. 



o 




100 200 300 400 500 
Time before present/Myr 

Figure 2. Time series models. The black line sliows the PNB 
model with T = 100 Myr and 13 = 0. The red line shows the QPM 
model with T = 100 Myr, /3 = 0, = 0.5 and Tq = 200 Myr; The 
blue line shows the SP model with A = 50 Myr and to = 200 Myr. 

RB/RNB: Random model in which a set of N times is 
drawn at random from a uniform distribution ex- 
tending over the range of the data. A Gaussian 
with standard deviation cr = 10 Myr is assigned to 
each of these, and then a constant B added be- 
fore normalizing. This is the RB model. The RNB 
("random no background") model is just the spe- 
cial case of i? = 0, which produces a model which 
is similar to our discrete data. In practice we fix 
N and B and calculate the evidence by averaging 
over a large number of realizations of the model. 
Specifically, when modelling the B5 and B18 data 
sets we fix S = 0, and = 5, 18 respectively. 

PB/PNB: Periodic model of period T and phase P 
(model PNB). There is no amplitude parameter be- 
cause the model is normalized over the time span of 
the data. Adding a background B to this simulates 
a periodic variation on top of a constant extinction 
probability (model PB). 

QPM: A quasi-periodic model in which the phase is a 
sinusoid with amplitude Aq , period Tq and phase 9 
(it becomes the same as the PNB model if Aq = 0). 

SP: A monotonically increasing or decreasing nonlinear 
trend in the extinction PDF using a sigmoidal func- 
tion characterized by the steepness of the slope, A 
and the center of the slope, to. In the limit that 
A becomes zero the model becomes a step function 
at to J and in the limit of very large A becomes the 
uniform model. 

SSP: Combination of SP and PNB. 

OM(P)/SOM(P): The orbital/semi-orbital model 
with/without spiral arms, defined in Section [4.51 

3.4. Numerical calculation of the evidence 

The integral in Eqn. [T] is a multidimensional integral 
over the parameter space, and cannot be calculated ana- 
lytically. As iu iBailer- Jones, (,2011ah . we estimate it using 
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Table 2 

The mathematical form of the time series models and their corresponding parameters. Time t increases into the past and Pu{t\B, M) is 

the unnormalizcd extinction probability density predicted by the model. 



model name 


Pu{t\e,M) 




parameters 


Uniform 


1 




none 


RNB/RB 


EiLiA/'(«;M"-'^)+ 


B 


a, N, B 


PNB/PB 


l/2{cos[27r(t/T + /3)] + 


1}+B 


T, 13, B 


QPM 


l/2{cos[27rt/T + Aq cos(2wt/TQ) + /3] + 1} 


T, /3, Aq, Tq 


SP 


[1 + g(t-to)/X]-l 




A, to 


SSP 


PNB+SP 




T, /3, A, to 


OM(P)/SOM(P) 


n{r^(t),v^{t)) 




r^{t = 0), v^{t = 0) 



Table 3 

Range of parameters adopted in the model prior parameter distributions. Except for OM(P)/SOM(P), a uniform prior for all parameters 
for all models is adopted which is constant inside the range shown, and zero outside. The prior PDF of parameters in the 
OM(P)/SOM(P) model is Gaussian and specified by the uncertainties in the initial conditions. 



model name 


range of prior 


Uniform 


None 


RNB/RB 


a = 10 Myr, N e {5, 18}, B G {0, y=j:} 


PNB/PB 


10 < T < 100, < /3 < 27r, B e [0, 1] 


QPM 


10 < T < 100, < ;3 < 27r, < Aq < 0.5, 200 < Tq < 500 


SP 


-100 < A < 100, 100 <to < 500 


SSP 


10 < T < 100, < /3 < 27r, -100 < A < 100, 100 < io < 500 


OM(P)/SOM(P) 


initial conditions (see Table [51 



a Monte Carlo method, by drawing model parameters at 
random from the prior distribution and calculating the 
likelihood at each. If the set of N parameter draws is 
denoted {9}, then Eqn. [T] can be approximated as the 
average likelihood 



PiD\M) ^ l-J2p(D\e,M) 



(8) 



In the following simulations we adopt iV = 10 000 unless 
noted otherwise. 

4. MODEL OF THE SOLAR ORBIT 

We now reconstruct the orbit of the Sun around the 
Galaxy over the past 550 Myr. This is done by integrat- 
ing the Sun's path back in time through a fixed gravita- 
tional potential, described in Section l47T] (The dynamics 
are reversible because only gravity acts; energy is not dis- 
sipated.) It has often been assumed that the solar orbit is 
periodic with respect to crossings of the Galactic plane 
and/or spiral arms; we investigate this numerically in 
Section [mi The stellar mass distribution corresponding 
to the potential gives the local stellar density which the 
Sun experiences in its orbit. In Section 14.51 we use this 
to derive the variation in the expected extinction rate. 

4.1. The Galactic potential 

We adopt an analytic Galaxy potential, ^g{R,z) (de- 
scribed in cylindrical coordinates) , comprising three com- 
ponents 

$G = *b + *h + *d • (9) 

The first two components are spherically symmetric dis- 
tributions which represent the b ulge and halo using 
Plmnmer's model ()Plummeiill911[ ) 



G M, 



h.h 



b.h 



(10) 



Table 4 

The para meters of Galactic pote ntial model (from 
IGarcfa-Sanchez et al.l (120011 ')) 

component parameter value 



Bulge 




= 1.3955 X 10^" 


Mq 




h = 


0.35 kpc 




Halo 


Mh -- 


= 6.9766 X 10" 


Mq 




bh = 


24.0 kpc 




Disk 


Md -- 


= 7.9080 X lO^o 


Mq 




ad = 


3.55 kpc 






bd = 


0.25 kpc 





and the third is an axi symmetric disk according to 
IMivamoto fc Nagail (|1975[) 



G Md 



+ {ad + ^Jz^^df 



(11) 



In Eqn. [10] and [TT] R is the galactocentric radius and 
z is the distance from the Galactic midplane. Mf,, 
and Md are the mass of bulge, halo and disk, respec- 
tively, ad and hd are the scale length and height (re- 
spectively) of the disk, and hd and hh are the scale 
lengths of the bulge and halo respectively. We adopt 
the numer ical value s for these as given in Table [4| (see 
(fGarcia-S anchez et al. 2001)). In addition to the these 
three components, we introduce into some of the mod- 
els a time-dependent potential due to the Galactic spiral 
arms, denoted ^s{R, (f>, z,t). It is defined below in Sec- 
tion [O 



There is, of course, significant uncertainty not only in 
the value of the parameters of the potential, but also in 
the functional form of our model. It is no doubt a sim- 
plification of the true potential of the Galaxy, so specific 
numerical values of quantities such as orbital periods and 
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amplitudes inferred should not be interpreted too liter- 
ally. Below we investigate the sensitivity of our model 
comparison to changes in the potential. 

4.2. Orbit calculation 

To calculate the motion of a body through the po- 
tential from given initial conditions, we solve Newton's 
equations of motion, which in cylindrical coordinates are 



9$ 

2RRcp^- — 
ocp 

oz 



(12) 



We solve these equations by numerical integration us- 
ing the Isoda method implemented in the R package 
deSolve, with a time step of 0.1 Myr. 

The initial conditions are the current phase space co- 
ordinates (three spatial and three velocity coordinates) 
of the Sun. These are derived from observations with 
a finite accuracy, so our initial conditions are Gaussian 
distributions, with mean equal to the estimated coor- 
dinate and standard deviation equal to its uncertainty 
(Table [5]). In order to calculate an orbit we draw the 
initial conditions at random from these prior distribu- 
tions, and a large number of draws gives us a sampling 
of orbits which will be used later (e.g. in the evidence 
calculations) . 

We derive our initial conditions from a number of 
sources: The distance to the Galactic centre comes from 
astrometric and spectroscopic observa tions of the stars 
near the black hole of the Galaxy (jEisenhauer et all 
I2003D . The Sun's displacement from the galactic plane is 
calculated fr om the photorn e tric ob servations of classical 
Cepheids by iMaiaess et all (|2009[ ). The Sun's velocity 
is calc ulated from Hipparcos data bv iDehnen fc BinnevI 

4.3. Spiral arms 

The model for the spiral arms is described by their 
geometry and their gravitational potential. However, for 
the arm crossing periodicity test in the next section we 
ignore mass of the spiral arms when calculating the solar 
orbit and only consider their location. Likewise, in one 
of the class of variants of the orbital model, OM and 
SOM (defined later), we ignore the arms entirely (for 
both the orbit and stellar density calculations). This is 
done so that we can see the additional effect of the arms, 
the form and mass of which are poorly determined by 
current observations. 

The geometric model comprises two logarithmic spi- 
ral arms, the positions of which in circular coordinates, 
(i?, ip), are given by 



0,(i?) = alog(i?/i?„ 



(13) 



where a is a winding constant, Rmin is the inner radius 
and (jjmin is the azimuth at that inner radius. The radius 
of the spiral arm ranges from Rmin to Rmax- Of th e 
various arm models offered bv iWainscoat et all ()1992D . 
we selected the main two spiral arms, 1 and 1', with (j)min 
given bv .VanhoUebcke et al., (,2009. ) and other parameters 



given by IWainscoat et al.l (|1992[ ) (see Table ^. Their 
location in the plane of the Galaxy is shown in the left 
panel of Figure |3l The arms rotate rigidly with constant 
angular velocity (pattern speed) of f ir, = 20kms~^kpc~^ 
(|Martos et al.ll200l lDrimme]|[2000[ ). 

The gravitational potential of the arms is described 
by the (first term of the) analytic potential of Eqn. 8 of 
ICox fc Gomea (l200l) . This is 



cos{N[cf>~MR,t)]) 



sech 



(14) 



where 



N 



R sin(a) ' 
l3i=KiH{l + 0AKiH), 
_ 1 + KiH + 0.3(gig)^ 
l + 0.3KiH 
cj)s{R,t)^c^,{R) + npt . 

The amplitude of the spiral density distribution is 

Pa{R, z) = pae-^sech^{z/H) , (15) 

where Rg is the radial length of the drop-off in density 
amplitude of the arms, po is the midplane arm density 
at fiducial radius ro , and H is the scale height of the spi- 
ral density. This modulated by a sinusoidal pattern in 
(f>, to give the overall density which corresponds approx- 
imateljO to the above potential 

Ps{R,(t>,z,t) ^ PAiR,z)cos[N{^- (f>,{R,t))], (16) 

where N is the number of arms and 4>s{R) is given in 
Eqn. 1131 The values of the relevant parameters are given 
in Table El 

4.4. Periodicity test 

In previous studies of the impact of astronomical phe- 
nomena on the terrestrial biosphere, it has frequently 
been assumed that the solar motion shows strict period- 
icities in its motion perpendicular to the Galactic plane, 
and sometimes also with respect to spiral arm crossings. 
We investigate this here using our numerical model. 

For each orbit k, we calculate the intervals between 
successive crossings, {At*}fc, (separately for midplane 
and spiral arm crossings), where i indexes the crossing. 
We then calculate the sample mean and sample standard 
deviation of these intervals for each orbit 



1 



Nk-1 



■i=i 



Nk-l 

E 

2 = 1 



(17) 



(18) 



The density corresponding strictly to the potential in Eqn. ll4l is 
rather messy, a nd is not necessary for this work. The approximate 
density in Eqn. 1161 is accura te enough for our purposes for small 
H/r and z/H not too large IICox fc U6me^[2002') . 
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Table 5 

The current phase space coordinates of the Sun, represented as Gaussian distributions, and used as the initial conditions in our orbital 

model 



R/kpc 


Vfl/kpc Myr~i 


(f>/va.d 


</>/rad Myr^i 


z/kpc 


Vz /kpc Myr-l 


mean 8.0 


-0.01 





0.0275 


0.026 


0.00717 


standard deviation 0.5 


0.00036 





0.003 


0.003 


0.00038 



Table 6 

The parameters of the geometric and potential model for the spiral arms. The parameters of the potential apply to both arms. 



geometric parameters 


arm 


a 


Rmin fkpc </)mi„/rad 


extent /kpc 


1 


4.25 


3.48 0.26 


6.0 


1' 


4.25 


3.48 3.40 


6.0 


potential parameters 


N 


Po/Mq kpc 


ro /kpc Rs /kpc 


H/kpc 


2 


2.5 X 10'^ 


8 7 


0.18 




Figure 3. The Solar orbit in the Galactic plane (left, thin lines) and perpendicular to the plane (right). The orbit in the left panel is in 
a reference frame rotating with the spiral arms (shown in thick lines). 



where Nk is the number of crossings in the fcth orbit. To 
assess the periodicity of the crossing intervals, we define 
the degree of aperiodicity as 

Ofc = CTfc/AtI . (19) 

An orbit with a = is strictly periodic. 

We investigated the variation of the aperiodicity of the 
solar orbit with the six parameters (initial conditions). 
This parameter space is too large to report on extensively 
here, but we find that the aperiodicity is most sensitive to 
R{t = 0) and (p{t = 0). In the following we vary these ini- 
tial conditions individually, by drawing 10'' samples from 
the corresponding initial condition distribution. (Larger 
sample sizes did not alter the results significantly) We 
simulate the solar orbits using the arm-free potential, 

4.4.1. Midplane crossings 

Some earlier studies claimed that Galactic midplane 
crossings trigger increases in terrestrial extinction due 
to an enhanced gamma ray or cosmic ray flux or due 
to larger perturbation of the Oort Cloud. These are 



directly related to the increased stellar density and in- 
creased occurrence of star-forming regions. The larger 
tidal forces are postulated to enhance the disruption of 
the P ort cloud (jRampino fc Stothersl[2000l : iMatese et all 
119951 ). and the higher density of massive stars — and 
thus high energy radiation as well as increased SN rate 
— raises the average flux the Earth is exposed to. The 
periodicity of the Sun's vertical motion - not least its 
period, phase and the assumed stability of this period - 
is central to these claims. We examine these using our 
model. 

The results of varying just the initial galactocentric 
radius of the Sun, R(t = 0), are shown in Figure S) We 
see in the top-right panel that about 90% orbits have an 
aperiodicity less than 0.1. In the lower two panels we see 
how a varies with the value of the initial condition and 
with the average crossing interval. 

The aperiodicity is 0.002 (nearly strict periodicity) at 
At = 61.8 Myr. This corresponds to a 1:1 resonance be- 
tween the vertical motion and the radial motion. Its 
value is close to a per iod in the biodi versity data of 
62 ± 3 Myr claimed by iRohde fc Minimi 12005.). Little 
should be read into this coincidence, however, as there 
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Figure 4. Periodicity test of midplane crossings varying just the current galactocentric radius of tlie Sun, R{t = 0). Top left: distribution 
of this initial condition in the simulations (it has a Gaussian distribution with parameters given in Table (SJl. Top right: cumulative 
probability of the aperiodicity parameter for the resulting orbits. Bottom left: the variation of aperiodicity with R{t = 0). Bottom right: 
the variation of aperiodicity with the average crossing interval, At^ . 



is no good (i.e. independent) reason to select the spe- 
cific initial condition that leads to this period over any 
other. Moreover, changing the parameters of the Galac- 
tic potential — which is not very well known — changes 
this period. (For example, if we increase the mass of 
the Galactic halo the values of At are decreased.) The 
other minimum in the aperiodicity in the bottom pan- 
els is 0.02 at At = 100.2 Myr. This corresponds to an 
approximately circular orbit in the midplane. If we set 
VR{t = 0) = 0, Vz{t = 0) = and z 0, this solar orbit 
would be strictly circular. 

The cumulative curve (top-right panel of Figure U) 
makes a sharp turn at a' = 0.1. This is because of 
a sudden decrease in the number of orbits with large 
aperiodicities. Similarly, the discontinuities in the lower 
panels are caused by changes in the (small) number of 
discrete plane crossings which occur for different aperi- 
odicity ranges. 

If we now vary the initial condition (/)(t = 0) instead, 
the periodicity test gives very similar results: we find 



a nearly strict periodicity at At — 60 Myr and another 
minimum in the aperiodicity at about 100 Myr. That 
means the nearly strict periodicity is mainly determined 
by a combination of R{t = 0) and (f>{t — 0). 

In summary, we see that the majority of the simulated 
orbits (90%) are quite close to periodic (a < 0.1) in their 
motion vertical to the midplane, although strict period- 
icity essentially never occurs. 

4.4.2. Spiral arm crossings 

Regarding spiral arms as regions of increased star for- 
mation activity and stellar density, the mechanisms of 
mass extinction considered for midplane crossing could 
likewise be applied to spiral arm crossings , and have been 
by so me authors (jLeitch fc Vasishtl 119981 : iGies &: Helsell 
120051 ). However, such studies have over simplified the 
solar motion by failing to take into account the consider- 
able uncertainties in the current phase space coordinates 
of the Sun and thus in its plausible orbits. Some stud- 
ies have even claimed a connection between spiral arm 
crossings and the terrestrial biosphere after having fit the 
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solar motion to the geological data, but such reasoning 
is clearly circular. 

We examine here the periodicity of spiral arm crossings 
(although we note that some studies in the literature 
claiming a spiral arm-extinction link just consider the 
crossing times and do not claim a periodic crossing). The 
crossing intervals are longer than with the midplane, so 
we only include in our analyses models in which there 
are at least three arm crossings. We assume that the 
arms have indefinite vertical extent, so that a crossing 
on the x-y plane is always a true encounter. In reality 
the Sun might pass over or under the arms, thus reducing 
the overall relevance of spiral arm crossings to terrestrial 
extinction. 

Figure [5] shows the result of this analysis for the 7407 
orbits (out of the original sample of 10 000) that exhibit 
at least three arm crossings. The cumulative probabil- 
ity (top-right panel) shows that about 40% of the orbits 
have an aperiodicity larger than 0.2. In other words, it is 
not very likely that the solar orbit and spiral arms are so 
tuned to give periodic crossings. The lower two graphs 
show how a varies with R{t = 0) and Atk- The numer- 
ous gaps in these plots are a consequence of the fact that 
not all orbits for certain ranges of R{t = 0) had at least 
three arm crossings, and so were removed from the anal- 
ysis. We see, therefore, that the crossing interval is very 
sensitive to R{t = 0). 

Note that we have neglected the mass of the arms in 
the orbital calculations. When we include it the values 
of aperiodicity increase and there is an even less clear 
dependence of a on i?(t = 0) or Atk- 

In summary, we find it unlikely that spiral arm cross- 
ings are even close to periodic. If the pattern speed of the 
spiral arms has not been constant in the past 550 Myr, 
or if the pattern itself has not been stable, then this con- 
clusion is strengthened further. 

4.5. Orbital model 

4.5.1. Derivation of the extinction rate from the stellar 
density variation 

As outlined in Section [1] various astronomical mech- 
anisms for biological extinction have been identified, 
including comet impacts (from Oort Cloud perturba- 
tion), gamma rays (from S Ne or GRBs), a n d cos - 
mic ra ys (from SN remnants; 'Ellis & Schramm' (1 19951 ): 
iGarcia-Sanchcz ct al. (2001); Gics & Helscl (2005)). The 
intensity of all of these depends on the local stellar den- 
sity. If we consider a general mechanism involving flux 
from nearby stars, then the flux from a single star is pro- 
portional to f /d^ , where / is the relevant surface flux and 
d the distance. The sum of this over the whole relevant 
volume of space around the Sun is proportional to the 
total intensity and thus the extinction probability (per 
unit time). 

Let us assume that the extinction rate, -E, is linearly 
proportional to the flux, and that the number density of 
relevant stars is proportional to the total stellar number 
density (stars per unit volume), n. Because the den- 
sity of spiral arms is much less than the density of the 
other components, we consider at flrst only the time- 
independent density arising from halo, disk and bulge. 
The density is calculated from the corresponding poten- 
tial (defined in Section |4T|) using Poisson's equation. 



In an axisymmetric cylindrical coordinate system, the 
extinction rate at the Sun is then 

E{Rq,Zq) = C I I "^•^.'^^ RdRdz 



d2 



-C 



n{R, z) 



{R-RqY + {z-zq) 



2 RdRdz, 



(20) 



where R and z are the galactocentric radius and height 
above the midplane, respectively, for some star, Rq and 
zq are the corresponding (time-varying) coordinates of 
the Sun, and C is a constant. Note that the stellar num- 
ber density, n, is proportional to the corresponding stel- 
lar density, p. Defining the distance from a star to the 
Sun as r = i? — Rq and Z = z — Zq, the extinction rate 
is 



E{Rq,zq)^C 



i{Rq +r,zo + Z) 
r^ + Z-^ 



{Rq + r)drdZ . 

(21) 



The flux from a star falls off as 1/d^, but we can trun- 
cate this integral at some upper distance because at some 
point the flux is too weak to influence the terrestrial bio- 
sphere. We take dth = 50 pc as an upper limit0 This 
is much smaller than the scale length of the disk and 
comparable to the scale height of the disk (see Table S]) , 
so we can approximate n(i?0 + r,ZQ + Z){Rq -\- r) by 
n{RQ, Zq + Z)Rq. The integral then becomes 



/ 

Integrating over r gives 



th — tlth 



n{RQ,ZQ + Z) 
r-^ + Z"^ 



drdZ 



(22) 



E{Rq,zq)^2CR^ 



7, ^ arctan(rfth/^) 
n(i?Q, Z+Zq) dZ . 

(23) 



The geometric factor ^'''^^an^dth/.g) j^jj-Qps close to zero at 
about Z — 2b pc, and does so much more rapidly than the 
stellar density term, which follows the vertical profile of 
the disk (which has a much larger scale height of 250 pc). 
Thus to a reasonable degree of approximation we can set 
n(i?o, Z -\- Zq) ~ n(i?0, Zq) in this integral. The integral 
is then just over the geometric factor, which gives some 
constant (dependent on c?th, but of no further interest). 
Thus we are left with 



E{RQ,ZQ)c^C'RQn{RQ,ZQ), 



(24) 



for some constant C . For the solar motion, the relative 
variation of Rq is less than that of n{RQ,ZQ), so we have 



E{Rq,zq) oc n(i?0,Z0) 



(25) 



In other words, the extinction rate is just proportional 
to the stellar density at the location of the Sun. The 



^ In the case of SNe. lEllis^& : SchramrJ II1995I ) conclude that only 
those which come within 10 pc of the Sun would have a significant 
impact on terrestrial life. GIlBs up to 1 kpc or even more could still 
have an effect on the Earth, but we ignore these because the GRB 
rate ( at low redshifts) is comparatively low (e.g. IDomainko et afl 

mm)- 
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Figure 5. Same as Figure|4] but now for the spiral arm crossings. 

approximations in Eqns. I22H25I still hold when we in- 
clude the low density spiral arms defined in Section |4.3[ 
in which case we must also introduce the explicit depen- 
dence on azimuth and time 

E{RQ,(f>Q,ZQ,t) (X n{RQ,(j)Q,ZQ,t) . (26) 

In the above model we assumed that the extinction 
rate is proportional to d~^, i.e. the influence falls off like 
a flux on the surface of a sphere. We could generalize 
this dependence to be d"*^/^ for fc > in order to reflect 
other mechanisms, e.g. tidal effects. 

In order to test the validity of the above approxima- 
tions, we compare in Figure|6]the extinction rate as given 
by Eqn. [50] (by numerical integration) with the stellar 
number density n{RQ, zq). We plot over ranges of Rq 
from 5kpc to lOkpc and zq from — O.Skpc to 0.5 kpc, 
in accordance with the ranges covered by the simulated 
solar orbits. We normalize the extinction rate (and the 
stellar density) by setting its integral over Rq and zq to 
be unity. In the upper row of Figure [51 the difference be- 
tween the stellar density and the extinction rate reaches 
a maximum in the midplane (z© = 0); this is on account 
of the relatively large density gradient at Zq = 0. The 



maximum difference is only about 10% of the peak value 
of stellar density for all values of fc. In the lower row, the 
largest difference is at the lower limit of Rq . Note that 
the value of fc has very little impact. 

In practice, most of the simulated orbits spend most 
of their time in the region 7 < i?0/kpc < 9 and — 0.3 < 
Z0/kpc < 0.3, where the differences between local stellar 
density and extinction rate variation are even smaller. 
Thus to within a few percent, the stellar density at the 
Sun is a good predictor of the extinction rate. The time 
variation of this density is the time series model forms 
the basis for what we refer to as the "orbital models", 
the forms of which we now define. 

4.5.2. Definition of OM(P) and SOM(P) 

The orbital model "OM" is the orbital model which 
does not include the spiral arm at all, neither in the 
gravitational potential (for calculation of the orbits) nor 
in the stellar density (for the extinction rate calculation). 
The orbital model OMP does include the spiral arm in 
both senses. Thus both OM and OMP are internally 
self-consistent. 

Once normalized, E{t) is just the quantity P{t\9^M) 
in Section 13.31 (and it is normalized to give unit inte- 
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Figure 6. Comparison of the extinction rate calculated numerically with the stellar density at the position of the Sun. The top row shows 
the variation as a function of zq with Rq fixed to 8kpc. The bottom row show the variation as a function of Rq with zq fixed to 26 pc. 
The columns from left to right are for A; = 0, 2, 4 in the model for the dependence of extinction rate with distance. 



gral over the span of the data) . The parameters of OM 
and OMP are the initial conditions of the orbit, and the 
corresponding priors are the Gaussian distributions sum- 
marized in Table [SI Thus one orbit calculated from one 
draw of the initial conditions allows us to calculate one 
likelihood for these models (for given data set). Repeat- 
ing this and averaging the resulting likelihoods gives the 
evidence for that orbital model (see Section 15^ . 

For both of these models we consider four variations, 
labeled 1-4, according to which initial conditions we vary 
(and therefore sample over to build up the set of orbits). 

In addition to these models, we define the "semi-orbital 
model" , SOM. This is derived from the OM simply by 
subtracting from the predicted extinction rate a constant 
value, ft., and setting all resulting negative values to zero. 
Here we simply set h to be the minimum value of the 
extinction rate (see Figure [7]). This is intended to model 
the situation in which the flux causing the extinction 
must rise above some threshold before it has an effect. 
(We might consider this as an adaption of life to the ex- 
traterrestrial flux background.) In analogy to OM, SOM 
excludes the spiral arm. SOMP is SOM with the spiral 
arm potential and density included. Once again we will 
consider four varieties according to which initial condi- 
tions are varied. 

5. RESULTS 

5.1. Evidences 

We now calculate the Bayesian evidence (Eqn. [T|) for 
the various models for each data set. This is done by sam- 
pling from prior probability distributions of the model 
parameters (P(0|M), Table [3]), calculating the likelihood 



(Eqn. m for discrete time series, Eqn. [7] for continuous 
time series) and then averaging these for that model and 
data set Eqn. [H 

To calculate the evidences for RNB and PNB models 
for the B5 and B18 data sets, we adopt a Monte Carlo 
sample size of 10^. In all other cases we use a sample 
size of 10^. Larger sample sizes did not alter the esti- 
mated evidence significantlyQ This sample size is given 
according to the sensitivity test of the evidence to the 
sample size for all models and all data sets. This test 
shows that the evidence estimated from 10^ draws in the 
prior distribution is close to the real evidence when there 
is a background either in the model or in the data set. 

As the absolute value of the evidence is not of interest, 
we report the ratio of evidence, the Bayes factor. Here we 
report Bayes factors with respect to the Uniform model. 
We regard a model as being significantly better than an- 
other when i ts evidence exceeds that of the oth er by a 
factor of 10 (|Jeffrevslll96ll : IKass fc Raftervl[T995h . Note 
that it is only meaningful to compare evidences — and 
therefore Bayes factors — for a fixed data set. 

The results are shown in Table [71 For the reference 
models we evaluate the evidence by sampling over all 
their model parameters, but in the case of OM and SOM 
we sample over just some of the parameters (initial con- 
ditions), keeping the others fixed, in order to investigate 
the impact of the different parameters. As shown in Sec- 
tion [44l the periodicity of solar orbit is most sensitive to 

^ For all the data sets, the standard error of the Monte Carlo 
estimates of the evidence is < 1% for OM models, < 3% for SOM 
and other models, and < 25% for RNB and (S)OMP models. 
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Figure 7. Orbital model (OM) and semi-orbital model (SOM). The left panel shows the variation of the local stellar density - and thus 
the extinction probability per unit time — in one particular orbit calculated from OM. The horizontal dashed line is a threshold, h, for 
truncating the stellar density to a minimum level, which gives rise to the extinction probability per unit time plotted in the right panel. 



Table 7 

Bayes factors and maximum likelihood ratios of the various time series models (rows) relative to the Uniform model for the various data 

sets (columns). OM(P)l-4 refer to the OM(P) model in which different initial conditions are varied: R{t = 0), 4'{t = 0), 
{R{t = 0), 4>{t = 0)}, {R{t = 0), VR(t = 0), 4>{t = 0), Vz{t = 0)}, respectively, and likewise for the SOM(P)l-4 models. The other initial 
conditions axe kept fixed. The RB and RNB models are intrinsically discrete, so are not applied to the two continuous data sets. 
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the initial conditions R(t = 0) and <f){t ~ 0). Wo there- 
fore calculate the evidence for the OM (and SOM) models 
with four different sets of initial conditions being varied: 
R{t = 0) only; (j){t = 0) only; {R{t = 0) and ^{t = 0)}; 
{R{t = 0), Vnit = 0), <j){t = 0), and = 0)}. In ah 
cases we fix cf){t = 0) and z{t = 0), the former because it 
has no impact on the solar motion in this axisymmetric 
potential, and the latter because the uncertainty in the 
current z position of the Sun has a limited impact on the 
subsequent orbit. To assess the effect of the spiral arm 
perturbation on the Bayes factors, we have selected four 
perturbed orbital models, OMPl, 0MP4, SOMPl and 
S0MP4, to compare with corresponding unperturbed or- 
bital models. 



For the B5 data set, the Bayes factors of all time sc- 
ries models relative to the Uniform model arc less than 
10. Thus none of these models arc a significantly better 
explanation of the data. One model, RNB, has a Bayes 
factor less than 0.1, indicating that we can discount this 
one as being an unlikely explanation. Given that the 
Uniform model is the simplest model of the sot, the prin- 
ciple of parsimony suggests that wc; should be satisfied 
with it as explanation. This does not deny the possi- 
bility that some other model shows significantly higher 
evidence. After all, we can only ever make claims about 
models which we explicitly test. 

The B18 data set includes more extinction events than 
the B5 data set, and not surprisingly it discriminates 
more between the models (the Bayes factors show a larger 
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Figure 8. Results for the B18 data sets. The upper panel shows 
the log (base 10) Bayes factor of the various models relative to 
the Uniform model. The lower panel shows the log (base 10) of the 
maximum likelihood ratio of various models relative to the Uniform 
model. 

spread) . (These results are also shown graphically in the 
upper panel of Figure [H) The OM models are favored 
somewhat more than the other models - e.g. the Bayes 
factor of 0M3 to SP is 0.63/0.019 = 33 - although again 
no model is favored significantly more than the Uniform 
model. In contrast, several models are significantly disfa- 
vored (RNB, SP, OMPl, 0MP4, SOMl, S0M2, SOMPl 
and S0MP4). In particular, the perturbed orbital mod- 
els, including OMPl, 0MP4, SOMPl, and S0MP4, are 
less favored by the data than their corresponding unper- 
turbed orbital models. All the other perturbed orbital 
models (not listed in Table [7]) also have lower Bayes fac- 
tors than the unperturbed orbital models. 

For the two continuous time series, RM and AOS, the 
difference between the evidences of all of the models is 
not significant^ Our broad conclusion is that no model 
significantly outperforms the Uniform model on any of 

As the RNB and RB models are obviously conceptually inap- 
propriate models for continuous data sets, we do not apply them to 
the AOS or RM data sets, so these values are missing from Table [7] 



the data sets. On the contrary, a few can be "rejected" on 
the ground of a significantly lower evidence. Recall that 
all of these models are predicting the extinction prob- 
ability (per unit time). In terms of the discrete data 
sets, the Uniform model just means that the mass ex- 
tinction events occur at random in time. We find this 
to be no less probable than a periodic or quasi-periodic 
variation of the probability, or a monotonic trend in the 
probability, etc. In terms of the continuous data sets, 
we obviously do not believe that the Uniform model is 
a good explanation of the clearly apparent variations in 
the extinction rate (see Figure [T]). But the analysis does 
tell us that this is no worse an explanation than the more 
complex models of the variation considered, such as peri- 
odic, orbital-model based etc. Clearly there must be yet 
other models which could explain the data even better. 
This may explain why previous authors have found an 
apparent periodicity in the data: the periodic model can 
explain the data to some degree, but actually no better 
than simpler models. 

5.2. Likelihood distribution 

We have seen that the evidence hardly discriminates 
between any of the models on the continuous data sets, 
and only between some of them on the discrete data sets. 
(This is by no means inevitable. In other problems the 
evidence can vary enormously between models.) This 
means that, on average over their parameter space, the 
models differ little in their predictions. It is nonethe- 
less interesting to see how the likelihood varies over the 
parameter space. (We would do this in particular to 
find the best fitting parameters, although these are only 
meaningful if the overall model has been identified as the 
best explanation of the data.) We focus here mainly on 
the PNB and OM models for the B18 data. We again 
normalize the likelihood for a model by dividing it by the 
likelihood of the Uniform model to form the likelihood 
ratio. As the latter model has no parameters, it is likeli- 
hood constant and equal to its evidence. The maximum 
value of the likelihood ratio we denote as MLR. 

Figure |9] shows how the likelihood varies over the two- 
dimensional space formed by the two parameters, pe- 
riod and phase, of the PNB model. There is significant 
variation. We see numerous local maxima, the largest 
likelihoods being around {T/Myr, /?/rad} — {60,4.5}. 
However, these maxima are rather narrow, so once the 
(much lower) likelihood in the other (equally plausible) 
regions are taken into account, the overall evidence for 
the model is not particularly high. If we are interested 
in the variation of likelihood with just the period, then 
we can marginalize this diagram over phase, and plot 
with respect to period, thereby forming a (Bayesian) pe- 
riodogram (lower panel). We see a clear peak around 
60 Myr. Thi s is coincident with the period of 62 ± 3 Myr 
identified bv lRohde fc Mulled (|2005D . It is tempting (but 
incorrect) to associate this peak value of the likelihood 
with the periodic model as a whole, and use it to claim 
a larger evidence for the periodic model. Certainly there 
is a degree of arbitrariness in the prior parameter dis- 
tribution — in this case a uniform distribution — and 
narrowing this range around this peak would clearly in- 
crease the evidence. For example, if we truncate the 
period range from its current value of [10, 100] Myr to 
[50, 80] Myr, then the Bayes factor relative to the Uni- 
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Figure 9. Likelihood distribution for the PNB model on the B18 
data set as a function of period and phase (upper panel) and period 
only (lower panel). In both panels we show the likelihood ratio of 
the PNB to the Uniform model, on the left as log (base 10) on a 
color scale. 



form model increases from 0.62 to 1.5. This is a rather 
modest increase, but we could increase it to a signifi- 
cantly high value with an even narrower prior. However, 
we may not use the data to find the best fitting param- 
eters and then claim that we should only consider the 
model near to these. We would need some other rea- 
soning or independent data for making such a selection. 
(The Rohde & Muller (2005) time series is not indepen- 
dent of B18 because both are based on the same paleon- 
tological data.) We do not see how, a priori, we could 
limit the plausible periods of periodicity to something as 
narrow as 50-80 Myr, let alone the much narrower range 
required to favor PNB significantly over other models. In 
the extreme limit of an infinitesimal region around the 
maximum likelihood, we end up doing model compari- 
son using the maximum likelihood. Just out of interest, 
these values are shown in Tableland plotted in Figure[51 
If we were to use this (incorrect) metric, then PNB and 
some of its variants have significantly higher likelihood 
than the Uniform model and several of the other models 



(although barely more than a factor of ten above the ran- 
dom model, RB). Another way of seeing why this is the 
wrong approach was already discussed in Section [3T] by 
focusing on the best fits we simply favor the more com- 
plex model. We could always define a more flexible model 
and so fit even better. 

We labor this point because many of the claims for 
a periodicity in biodiversity data have made use of a 
maximum likelihood approach (of which is a special 
case) or something equivalent. We must instead use the 
evidence for model comparison. (Maximum likelihood 
may be used for estimating the best parameters once we 
have established we have the best model.) If the peri- 
odic model were in fact the true one, then of course only 
one period and phase would be true. In that case the 
likelihood around these values would be so high as to 
result in a large evidence even when averaging over the 
br oader parameter spa ce (see simulations in Section 4.1 
of iBailer-JonesI ()2011af ) for a demonstration). 

Incidentally, the fact that we find a dominant period 
at all in the B18 data set is actually not that unlikely. 
The (Bayesian) periodogram of samples drawn at ran- 
dom from a uniform distribution often exhibits a period 
which has a likelihood larger t han that of the true Uni- 
form model (see Section 4.2 of iBailer- Jones! (|2011al) ). In 
other words, it is often possible to explain a random data 
set with some period, which is just a testament to how 
flexible the periodic model is. 

Moving on from the periodic model, we show in Fig- 
ure [To] the likelihood distribution for OMl and 0M2, 
i.e. where we vary the initial conditions R{t = 0) and 
4){t = 0), respectively. The likelihood ratio varies by a 
factor of up to lO'*, but its absolute value is never more 
than about two. That is, no value of the intial condi- 
tions gives a model much more favorable than the Uni- 
form model, whereas as some are far less favorable. If we 
had lower uncertainties on these phase space coordinates 
of the Sun, then we might be able to conclude something 
more definitive. For example, if R(t = 0) = 7.5 kpc then 
the OM model would be even less favored. We see a sim- 
ilar degree of variability for the other OM models and 
data sets listed in Table [T] 

We performed a similar analysis for the other model, 
but for the sake of space report only the maximum like- 
lihood ratios in Table [7l 

For the B5 data set, the RNB model has the highest 
maximum likelihood ratio, yet its evidence was the low- 
est. This indicates that while one particular instance of 
RNB fits the data well, overall it is a poor model. 

For the RM and AOS data sets, the evidences are very 
similar for all models. This means that the data are 
not able to discriminate between these models very well: 
they are equally good (or bad). However, the time se- 
ries analysis model used here is not best suited to these 
data sets. These data can be better interpreted as valued 
time series, ones in which we have an extinction magni- 
tude attached to each time (both, in general, with un- 
certainties) , rather than a tim e variable probability of 
extinction. IBailer- Jones! (|2012| ) has extended the present 
model in order to work with such data sets; the results 
of its application to RM and A08 will be reported in a 
future publication. 

In summary, we find that for none of the data sets is 
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Figure 10. The variation of initial conditions (left column) and corresponding variation of likelihood ratio relative to the Uniform model 
(right column) for the OM model on the B18 data set. The top row is for OMl, in which only R{t = 0) is varied. The bottom row is for 
OM2, in which only if>{t = 0) is varied. 



any model particularly favored over the simple uniform 
one. In particular, there is no evidence to suggest that 
the orbital model (with or without a background extinc- 
tion level) is a particularly good explanation of the data. 

5.3. Sensitivity test 

The evidence is of course sensitive to the choice of prior 
parameter distribution, and often we have little reason to 
make a very specific choice. Here we test the sensitivity 
of the evidence to this, as well as to the parameters of 
the Galactic potential used in the orbital models, and to 
the age uncertainties for the discrete data sets. 

The age uncertainties in the discrete data are taken 
into account by the likelihood function. However, it is 
often difficult to estimate uncertainties, and we addition- 
ally made a plausible, but not unique, translation of the 
estimated duration of a stratigraphic substage in order 
to estimate uncertainty (which is the standard deviation 
of a Gaussian for each event; see Section [2.1^ . To see 
how this affects our results, we scale the age uncertain- 
ties in the B18 data set by a constant factor of 1/4, 1/2, 
2, and 4. For each of these modified data sets we cal- 



culate the evidences for the models (S)0Ml-4 and the 
Uniform model and recalculate the Bayes factors rela- 
tive to the Uniform model. These are plotted in the top 
panel of Figure [TT] The Bayes factors change by just a 
few percent, so a precise age uncertainty is not necessary. 

As a second test, we scale in the same way the uncer- 
tainties of the initial conditions of the orbital models (i.e. 
we change the width of the prior parameter distribution) . 
The results are shown in the bottom panel of Figure 1111 
The change in evidence for any particular model is larger 
than in the previous case, but in most cases less than a 
factor of 5 (except for the models including the spiral 
arms). Moreover, the absolute value of the Bayes factor 
remains below one. 

Our results are therefore also insensitive to consider- 
able imprecision in the uncertainties in the phase space 
coordinates of the Sun. This (and the previous conclu- 
sion) is also true for the B5 data set. 

As a third sensitivity test, we allow the number of sim- 
ulated random events, N, and the standard deviation of 
each event, a, in the RNB and RB models to vary. (Ear- 
lier we fixed N = 5 when drawing models for the B5 
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Figure 11. Test of the sensitivity of the evidence to variations in 
the data age uncertainties and model priors. Both panels plot (on 
the vertical axles) the values of the the Bayes factor of the varied 
OM(P) and SOM(P) models (relative to the Uniform model) for 
the B18 data sets. In the upper panel, the labels of the horizontal 
axis (Xsd) denote the varied B18 data sets, where Xsd means that 
we have multiplied the age uncertainties of the B18 time series by 
X for X=l/4, 1/2, 1, 2, 4. In the lower panel, the labels of the 
horizontal axis (Xsd) denote variations of the uncertainties in (i.e. 
spread of) the initial conditions, where Xsd means multiplying the 
uncertainties of the initial conditions by a factor of X (X=l/4, 1/2, 
1, 2, 4). For each value of Xsd wc have, for clarity, slightly offset in 
the horizontal direction the values for the OM(P) models (square 
points) from those for the SOM(P) models (diamond points). 

data set and iV 18 for the B18 data set.) For the B18 
data set, we find that a larger number of peaks or larger 
standard deviation in the RNB model produces a signif- 
icantly larger Bayes factor (see Table although it is 
still below unity. The RB model shows much less sensi- 
tivity to N and a. We similarly recalculate the evidence 
for the other models in response to various perturbations 
of their priors, as also listed in Table [51 The resulting 
Bayes factors do not change by more than a factor of 10 
in any case, and often by much less. 

Finally, we test the sensitivity to the Bayes factors to 
changes in the parameters of the Galaxy model (canon- 
ical values listed in Table |4|) . The results are shown in 



Table 8 

The Bayes factors (relative to the Uniform model) on the B18 
data set for models with priors varied. Each prior is varied 
individually (listed in the middle column) with the other fixed at 
their canonical values. 



models 



PNB 



PB 



QPM 



RNB 



RB 



SP 



SSP 



varied prior 



none 
50 < r < 80 
10 < T < 200 
10 < T < 400 



none 
B = 1/2 
B = 2 



none 
0<Aq< 1/4 

0<Aq <1 
100 <Tq < 300 
100 <Tq < 500 



none 
r = 5 Myr 
= 20 Myr 
N = 9 
N = 36 



none 
= 5 Myr 
= 20 Myr 



B = 
B -. 



/27r<T 
2 



✓27rCT 

N=9 
N=36 



none 
-200 < A < 
< A < 200 
10 < to < 500 



none 
10 < T < 200 
-200 < A < 200 
10 < to < 500 



BF 



0.62 
1.5 
0.31 
0.15 



0.80 
0.88 
0.97 



0.85 
0.62 

0.54 
0.61 
0.58 



0.00050 
3.7 X 10~" 
0.026 
0.0085 
0.29 



0.40 
0.73 
0.25 
0.13 
0.45 
0.84 
0.20 



0.019 
0.070 
0.10 
0.047 



0.18 
0.17 
0.37 
0.27 



Table |9|) for the B18 data set. If we double the mass 
of the halo, for example, then the evidence for the OM 
and SOM models changes by no more than a factor of 
three. Some other changes produce smaller effects, some 
larger, but not more than by a factor of 5 (and note that 
a change in a factor of two of the scale lengths is be- 
yond what is consistent with observed data) . Changes in 
the parameters of models which include spiral arms (the 
OMP and SOMP models) can produce larger changes in 
the Bayes factors. However, most significantly, none of 
these changes produce a Bayes factor greater than one. 
In other words, none of these changes result in the orbital 
model becoming a better explanation for the paleonto- 
logical than the Uniform model. 

In summary, we find that the evidences for most mod- 
els are not particularly sensitive to the age uncertainties. 
Galaxy model parameters, or reasonable changes to the 
prior parameter distributions. 

5.4. Testing the discriminative power 

Here we investigate how well our analysis can discrim- 
inate between models, by using simulated data which 
have been drawn from one of these models. Given suffi- 
cient data, such a discrimination will always be possible 
to some threshold Bayes factor, but here we are inter- 
ested in the case where the data have similar properties 
(in particular, sparsity) to the real data we have been 
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Table 9 

The Baycs factors on the B18 data set for orbital models with varied Galaxy parameters. 



variation 


UMi 




UMo 




UMr 1 


UMr4 


bUMi 


0/^T\ /TO 


0/^T\ /TO 


bUM4 


SOMPl 


S0MP4 


none 


0.74 


0.72 


0.63 


0.65 


0.014 


0.022 


0.051 


0.037 


0.032 


0.032 


1.3x 10~ 


-4 


1.2x 10" 


-3 


2Mb 


0.58 


0.45 


0.54 


0.55 


0.011 


0.018 


0.043 


0.011 


0.027 


0.028 


5.8x 10" 


-4 


3.8x 10" 


-4 


1 /2Mb 


0.78 


0.85 


0.67 


0.67 


0.0066 


0.021 


0.040 


0.045 


0.030 


0.030 


1.4x10" 


-4 


1.9x 10" 


-3 


2bt 


0.74 


0.72 


0.64 


0.63 


0.0036 


0.013 


0.055 


0.042 


0.033 


0.032 


4.6x 10" 


-5 


6. Ox 10" 


-4 


l/2bb 


0.74 


0.73 


0.65 


0.63 


0.019 


0.025 


0.051 


0.037 


0.032 


0.030 


2. Ox 10" 


-4 


1.5x 10" 


-3 


2Mh 


0.60 


0.68 


0.57 


0.57 


n r\r\r\-\ a 

0.00014 


0.0031 


0.080 


0.069 


0.087 


0.083 


4.D X lU 


6 


y.4x lu 


-5 




U.ol 


U.Do 


U.DD 


U.DO 


U.Ull 


U.U14 


n on 
U.zU 


u.uoy 


U.lo 


U.lo 


2.3x10" 


-4 


2.8x10" 


-4 


2bh 


0.81 


0.52 


0.66 


0.65 


0.011 


0.0061 


0.19 


0.027 


0.15 


0.15 


1.1x10" 


'3 


3.5x10" 


-A 


l/2bh 


0.073 


0.048 


0.11 


0.11 


0.024 


0.062 


0.0046 


0.0012 


0.010 


0.010 


7.5x10" 


-3 


1.4x10" 


-2 


2Md 


0.21 


0.063 


0.33 


0.32 


0.035 


0.058 


0.0047 


0.0046 


0.027 


0.024 


2.2x10" 


'5 


7.6x10" 


-A 


l/2Md 


0.59 


0.56 


0.49 


0.48 


0.0031 


0.0037 


0.0073 


0.0018 


0.0086 


0.0088 


3.3x10" 


-4 


9.6x10" 


-A 


2ad 


0.86 


0.89 


0.76 


0.76 


0.13 


0.12 


0.10 


0.088 


0.069 


0.067 


6.3x10" 


-3 


7.0x10" 


-3 


l/2ad 


0.63 


0.66 


0.55 


0.55 


0.021 


0.019 


0.015 


0.011 


0.022 


0.026 


2.5x10" 


'3 


9.1x10" 


-4 


26d 


1.1 


1.1 


0.93 


0.93 


0.0056 


0.020 


0.028 


0.030 


0.027 


0.025 


1.3x10" 


-A 


4.2x10" 


-3 


l/26d 


0.72 


0.87 


0.60 


0.56 


0.00085 


0.0073 


0.15 


0.17 


0.10 


0.090 


1.3x10" 


-A 


1.5x10" 


-3 



using. 

We investigate this by simulating a number of of time 
series from tlie RNB, OMl, and PNB models (which we 
will refer to as "generative models" when used in this 
way, to distinguish from their use to calculate the evi- 
dence on given data). For each generative model, we fix 
the parameters to certain values, then sample 18 events 
from the resulting P{tj\6,M) to give a simulated time 
series, to which we then attach the measured age uncer- 
tainties. We generate ten time series in this way (and 
below we average the Bayes factors over these and re- 
port that). For the OMl and PNB generative models we 
repeat this at ten different values of the solar initial ra- 
dius parameter (OMl) or period parameter (PNB) in the 
generative model. (RNB has no parameters). We repeat 
the whole process a second time but using simulated age 
uncertainties drawn from a log normal distribution with 
standard deviation and mean calculated from the mea- 
sured age uncertainties. We finally repeat the process a 
third time for the OMl and PNB generative model, but 
now drawing data to have the same time sampling as our 
continuous data sets (for which age uncertainties are not 
used; see section [3.2. 2p . (We do not do this for the RNB 
model as it predicts discrete events.) 

For each simulated data set we calculate Bayes fac- 
tors for the RNB, OMl, and PNB models relative to the 
Uniform model. For the data generated from the RNB 
model (with age uncertainties taken from the data), the 
Bayes factors for the three models are as follows: 0.18 
for RNB; 0.57 for OMl; 0.52 for PNB. (We get almost 
identical values when the age uncertainties were drawn 
at random). Thus no model - not even the true one - 
is favored over the Uniform model (although none is sig- 
nificantly rejected either). This is not that surprising, 
however, because with only 18 events, and with the evi- 
dence effectively averaging the predicted times of events 
from the RNB model over all time, the Uniform and RNB 
models end up with similar predictive power. This is un- 
avoidable, because with the RNB model we cannot decide 
in advance where the events are: we must average over 
all possibilities. 

The results for applying the models to the data gen- 
erated from the OMl and PNB models are shown in 



Fig. [121 where the horizontal axis shows how the Bayes 
factor varies with the one parameter which is varied in 
these generative models. The top row shows the results 
for data drawn from the OMl model, for the discrete 
(B18-like) data (left) and the continuous data (right). 
We see that the (true) OMl model is not significantly 
favored in either case (Bayes factor always less than 10), 
although not disfavored either (Bayes factor more than 
0.1). In particular, the continuous data show no discrim- 
inative power. 

The lower two rows show conceptually the same thing, 
but now for data drawn from the PNB model for two dif- 
ferent values of the phase parameter (the two rows). Here 
we see that, at least for longer periods, the PNB model 
is generally correctly identified (on the basis of a large 
Bayes factor), when using the discrete data sets. Yet the 
continuous data still show no discriminative power. 

This difference between discrete and continuous data 
sets is not unexpected. In the former, the likelihood is a 
product of the likelihood for many events, each of which 
is the convolution of the event with the model. In the 
latter, the likelihood is just the result of a single convolu- 
tion of a continuous model over continuous data. This is 
not the best approach for modeling continuous data. A 
better choice is the recently develo ped continuous tini e 
series modeling method described in lBailer-JonesI ()2012| 1 , 
which will be used in future work. 

Clearly one could perform many more tests with more 
simulated time series, varying different parameters in the 
generative models and with different permutations of the 
values of the fixed parameters. No doubt there are parts 
of parameter space where some models are favored over 
others, in particular if we adopted more informative pri- 
ors. Thus while these results on simulated data give some 
check on the discriminative power of the method and 
data, they should not be over-interpreted to say anything 
too general. Nonetheless, the tests we have done confirm 
what we concluded based on the analysis of the real data. 
Specifically, while our analysis of the real data does not 
allow us to claim evidence in favor of the orbital-based 
models, it also cannot rule out these models. This is due 
partly to the lack of predictive power of the data, and 
partly to the large fiexibility (or broad prior parameter 
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space) of the models. Better constraints on the solar 
orbit would help reduce the latter. 

6. SUMMARY AND CONCLUSIONS 

We have used a Bayesian model comparison method 
to examine how well different time series models explain 
the variation of biodiversity over the Phanerozoic eon 
(the past 550 Myr). One class of models is derived from 
the orbit of the Sun around the Galaxy, which we recon- 
structed from a model of the Galactic mass distribution. 
Our model comparison takes into account uncertainties 
in the data as well as uncertainties in the reconstructed 
path of the Sun. We have compared the evidence for this 
model with that of various other reference models of no 
particular causal origin. All models are stochastic in the 
sense that they predict only the time variation of the ex- 
tinction probability, rather than the exact magnitude of 
the extinction rate or the times of mass extinctions. 

As part of this analysis we investigated the properties 
of plausible solar orbits (i.e. those consistent with the 
accuracy of the present phase space coordinates of the 
Sun). We find that the majority of orbits have a motion 
perpendicular to the Galactic plane which is not far from 
periodic, although a precise period cannot be inferred 
due to the uncertainties in the present solar phase space 
coordinates as well as the exact mass distribution (gravi- 
tational potential) of the Galaxy. Thus any claims which 
try to link a variation in geological biodiversity, cratering 
or climate records to this solar motion must consider the 
motion as quasi-periodic rather than strictly periodic. 

In contrast, only about half of the simulated orbits 
showed periodic spiral arm crossings, even for a very sim- 
ple, rigidly rotation arm model. Indeed, many of the or- 
bits did not encounter the spiral arms more than once. It 
should be noted that the shape and pattern speed of the 
spiral arms is poorly known (and they may not even be 
long-lived) , so any claims of a causal connection between 
spiral arm passages and terrestrial conditions should be 
treated with due skepticism. 

We have shown how the evidence (marginal likelihood) 
should be used to do model comparison, as opposed to 
selecting the model which gives the best single fit. The 
reason is that an arbitrarily complex model can always 
be tuned to fit the data arbitrarily well, yet that does not 
make it a good model. By averaging the likelihood over 
the parameter space, the evidence uses the rules of prob- 
ability to trade off the quality of the fit with the model 
plausibility in a quantitative fashion. Of the models in- 
vestigated, we do not claim any one of them to be "true" . 
Indeed, no model is exactly true in reality. All we can 
hope to do is to find the best of the ones tested so far. 

We find that none of the models tested — including 
periodic, quasi-periodic and orbital-based — explain the 
discrete data sets better than a Uniform model. In other 
words, the time distribution of mass extinction events 
is consistent with being randomly distributed in time. 
There is no need to resort to anything more exotic. 

The Uniform model is also no worse than other mod- 
els for the continuous data sets. This does not mean 
that we believe the extinction rate has been constant 
over the Phanerozoic, but rather that none of the other 
(more complex) models is significantly better. Assuming 
the variations in extinction seen in Figure [T] are true — 
something we have no reason to doubt — then this tells 



us that there must be other models, not yet tested, which 
could explain the data better. This will be investigated 
in future work using a model more suited for these types 
of time series. 

We found in particular that the orbital-based extinc- 
tion model is not favored by the data. This conclusion 
is robust to changes in the parameters of the Galaxy 
model and to the magnitude of the uncertainties of both 
the solar phase space coordinates and the ages of the ex- 
tinction events. On the other hand, our analysis of simu- 
lated data showed that even if the orbital model were the 
true one, our analysis could not have identified it with 
either the discrete or (in particular) the continuous data 
sets. This ultimately comes down to a combination of 
a lack of discriminative power in the data, plus a large 
fiexibility (or prior parameter space) in the models. Of 
course, if the orbit of the Sun could be much better de- 
termined then it is possible that this model would then 
be more — or less — favored by our analysis. We remind 
the reader that our orbital model adopted an extinction 
mechanism in which the extinction rate is proportional 
to the integrated "fiux" (of a non-specified type) from 
nearby stars. A radical change in this mechanism would 
of course correspond to a quite different model, which 
could give different results. Thus we do not claim that 
the solar motion plays no part in terrestrial extinction, 
nor that astronomical mechanisms are irrelevant. 

Indeed, it is quite plausible that the biological extinc- 
tion rate has been affected by many factors, and so any 
attempt to connect them solely to the solar motion, or 
indeed to any simple analytic model, is doomed from the 
start. We have addressed this to some extent by includ- 
ing compound models and the semi-orbital model, but 
clearly one could do more. However, given the present 
uncertainties of the reconstructed solar orbit, it seems 
unlikely that one could draw a strong conclusion on the 
positive relevance of the solar orbit on the basis of cur- 
rent geological data. This, indeed, is the main conclusion 
of this work, plus the confirmation that periodic models 
are not a good (or necessary) explanation of the biodi- 
versity variation. There is some hope that, in the fu- 
ture, results fro r n the Gaia survey of the Galaxy (e.g., 
iLindegren et al.l ()2008D ) will improve our knowledge of 
the Galactic potential, spiral arms and inferred solar or- 
bit, to the extent that this study can be repeated to give 
conclusions of greater certainty. 
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